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We apply a biomass pyrolysis model, based on the model developed by Haseli et al. [4], which can be 
used in combination with Direct Numerical Simulation. The pyrolysis model is combined with a model for 
particle tracking to simulate 3D turbulent particle-laden channel flow with biomass particles undergoing 
pyrolysis in nitrogen. Transfer of momentum, heat and mass between gas and particles are fully taken into 
account. The effects of this transfer are analyzed and quantified in terms of the delay in the conversion 
or pyrolysis time. The delay is shown to depend on the initial volume fraction (number of particles) and 
on the size of the particles. The two-way coupling effects are relevant at volume fractions >10 -5 . For a 
fixed volume fraction, gas-particle interaction induces a delay in the devolatilization, decreasing with 
increasing particle size. Using this model, we also performed simulations of realistic biomass particle size 
distributions in order to compare two-way and one-way coupling. 

© 2014 Elsevier B.V. All rights reserved. 


1. Introduction 

Biomass has the attention of researchers and industrial compa¬ 
nies because it is believed to play an important role in the future 
since, unlike other renewable energies, it is available in large quan¬ 
tities in many regions of our planet [1], Secondly, but of equal 
importance, biomass research is important in view of increasing 
concerns about the environmental impact of most other meth¬ 
ods of energy conversion. For instance, the target of reducing CO2 
emissions from coal-fired power plants can be accomplished with 
co-firing at high biomass fraction and high oxygen concentration 
[2,3]. 

During the pyrolysis of biomass, virgin biomass is converted into 
char and volatiles are released. The volatile composition, concentra¬ 
tion and thermal properties depend on the fuel type, temperature, 
pressure, heating rate and reaction time [4], Various models of a sin¬ 
gle biomass particle undergoing pyrolysis are available in literature. 
Haseli has presented an overview of the available models in [5]. 
These models are able to capture the variation of important particle 
quantities during pyrolysis. However, the complexity of these mod¬ 
els requires a computational cost which is unfeasible when dealing 
with a large number of particles undergoing thermo-chemical con¬ 
version and interacting with the surrounding turbulent gas flow. In 
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some applications, like the design of combustors or gasifiers, only 
a few parameters are of interest; for example, the mean tempera¬ 
ture or surface temperature of the particles, the amount of volatiles 
released during pyrolysis and the ignition time or the conversion 
time of the particles. For simulating biomass pyrolysis with a large 
number of particles, a simplified biomass pyrolysis model is nec¬ 
essary. Here, we present a biomass pyrolysis model based on the 
model developed by Haseli et al. [5,6], which has the advantage 
to be suitable in combination with Direct Numerical Simulation 
(DNS). We combined the pyrolysis model with a model for particle 
tracking to simulate 3D turbulent particle-laden channel flow with 
particles consisting of biomass undergoing pyrolysis and interact¬ 
ing with nitrogen. The turbulent gas flow is solved in an Eulerian 
framework and the particles in a Lagrangian way. To the best of 
our knowledge, this is the first 3D Euler-Lagrange formulation of 
biomass pyrolysis. 

The presence of particles in a turbulent gas requires that the 
flow around each particle is solved. Despite the enormous progress 
in computing power, it is still impossible to simulate millions of 
particles interacting with a turbulent flow up to all details of the 
flow [7-9], Therefore, in case of turbulent flow with large num¬ 
bers of particles with sizes smaller than the Kolmogorov scale jj, 
it is common practice to adopt a point-particle approach to keep 
the computational cost at acceptable levels [10,11], This approach 
allows numerical simulations with millions of particles [12] and 
has been used to perform direct numerical simulations (DNS) as 
well as large-eddy simulations (LES) in many applications [13-16]. 
By employing Haseli’s pyrolysis model, we extend this approach to 
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Nomenclature 

otg thermal diffusivity of virgin biomass (m 2 /s) 
ac thermal diffusivity of char (m 2 /s) 

A h p specific heat of pyrolysis (J/kg) 

At time step (s) 

q heat flux applied to particle surface (W/m 2 ) 

€ emissivity of particle 

li g dynamic viscosity of gas (kg/(m s)) 
v kinematic viscosity of gas (m 2 /s) 

<f> particle volume fraction 

Re r frictional Reynolds number 

p p particle mass density (kg/m 3 ) 

a Stefan-Boltzmann constant (W/(m 2 1< 4 )) 

x p particle relaxation time (s) 

r w wall shear stress (kg/(m s 2 )) 

u velocity (m/s) 

V, particle velocity (m/s) 

Xj particle position (m) 

C c specific heat capacity of char (J/(kg K)) 

C vo ; specific heat capacity of volatiles (J/(kg K)) 
dp particle diameter (m) 

h heat transfer coefficient (W/(m 2 K)) 

k g thermal conductivity of gas (W/(m K)) 

rrij particle mass (kg) 

R particle radius (m) 

r radial coordinate (m) 

r c char front (m) 

r t thermal front (m) 

T temperature (K) 

t time (s) 

To initial particle temperature (K) 

T p pyrolysis temperature (K) 

T s particle surface temperature (K) 

T w wall temperature (K) 

u T friction velocity (m/s) 

v c velocity of char front (m/s) 

Re p particle Reynolds number 

H half the channel height (m) 


biomass pyrolysis including the heat and mass exchange between 
gas and particle. Usually, this approach implies uniform temper¬ 
ature of the particle, whereas in the model presented here, the 
thermal diffusion inside the particle is taken into account although 
the particle is modeled as a point. 

In this paper the steps necessary to obtain the ordinary differen¬ 
tial equations from the work of Haseli et al. [6] and to describe the 
pyrolysis coupled to the flow and the subsequent steps to enable the 
coupling of the model to DNS of turbulent gas flow are described. 
The flow regime in burners is typically turbulent to enhance mix¬ 
ing and, therefore, the combustion efficiency. We test the biomass 
pyrolysis model in a turbulent channel flow. A real burner geom¬ 
etry was used by Ghenai and Janajreh [17], who performed a 2D 
simulation of co-firing biomass with coal. Ghenai and Janajreh 
[17] did not solve the turbulent fluctuations directly but used the 
RANS equations adopting the k-e turbulence model and a stochastic 
tracking model to predict particle dispersion. In our formulation the 
particle-gas interactions are modeled in detail. We restrict to the 
pyrolysis phase, whose prediction, in case of biomass, is very impor¬ 
tant for a good prediction of the subsequent combustion phase. 
Implementation of combustion will be a future extension of this 
research. 

An important parameter in the pyrolysis of biomass is the con¬ 
version time, i.e., the time biomass particles need to fully convert to 


char. Due to the very high content of volatiles in biomass, the con¬ 
centration of volatiles influences the combustion process. With our 
model we can compute the volatiles concentration at any time and 
take into account the effects of particle-gas interactions as well. We 
focus on the pyrolysis conversion or volatilization and in particular 
analyze the effect of the two-way coupling on the volatilization. 
We also simulate the mass loss of a realistic biomass particle 
size distribution to compare the differences between a two-way 
and a one-way coupling formulation and investigate under which 
conditions gas-particle interaction considerably delays biomass 
pyrolysis. 

The structure of the paper is as follows. In the next section 
Haseli’s model is introduced and the steps to make the model suit¬ 
able for implementation in DNS are described. In Section 3, the gas 
model is introduced. Then, the numerical methods used to perform 
the simulations and the simulation set up are described in Sections 
4 and 5, respectively. Section 6 is dedicated to the presentation and 
discussion of the results. 


2. Biomass particle model 

In this section, we present the pyrolysis model used to derive 
the equations for the particle temperature during pyrolysis. Then, 
the equation of motion used for tracking each single particle in the 
system is presented. 

The biomass type considered in this work is torrefied wood. Tor- 
refaction or drying is a thermo chemical treatment of biomass at 
200 to 320 °C. During this process, carried out in the absence of oxy¬ 
gen, the water contained in the biomass is released obtaining dry 
biomass. The torrefaction helps to produce much better fuel quality 
for combustion and gasification applications. It improves the grind- 
ability and the hydrophobicity, very important to store the biomass 
outside. The energy required for torrefaction is partly compensated 
by the energy reduction for grinding torrefied wood compared to 
raw wood. The energy requirement for grinding biomass chips is 
about 1/6 of the heating value. The torrefaction can reduce this 
energy up to 90% [18], 

After milling, two of the main differences between wood and 
torrefied wood are the smaller size and the more spherical shape 
of the torrefied biomass particle [19]. The biomass particle model 
presented here is based on the assumption of spherical shape of the 
particle, which is a better assumption for torrefied wood than for 
the raw material. 

2.1. The mathematical model of biomass particle pyrolysis 

The biomass pyrolysis model of Haseli et al. [5,6] consists of sev¬ 
eral phases for modeling the overall process. The pyrolysis process 
is only one of these phases, which are described in the following. 

Consider a spherical particle with a uniform initial tempera¬ 
ture. When a positive heat flux ij is applied to its surface, the 
surface temperature starts to rise and, due to the thermal con¬ 
ductivity of the particle, the temperature profile inside the particle 
starts to change. Theoretically, the velocity at which the thermal 
disturbance propagates through the particle is infinite. However, 
following Haseli et al. [5,6], immediately after the heat flux is 
applied, the particle is divided into two regions, the outer region 
closer to the surface, in which the temperature has changed, and 
the inner region around the center of the particle, in which the 
temperature is still unchanged (Fig. 1(a)). The surface which sepa¬ 
rates these two regions is called thermal front and moves from the 
particle surface towards the particle center. Although the thermal 
front is assumed to move with a finite velocity, the model has been 
demonstrated to predict results which agree well with a full ther¬ 
mal model. This method of introducing a thermal moving front is 
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Fig. 1. Stages of a biomass particle undergoing pyrolysis: heating of virgin biomass 
(a), pre-pyrolysis (b), pyrolysis (c), and post-pyrolysis heating (d). T 0 , initial particle 
temperature; T p , pyrolysis temperature; T s , particle surface temperature; r, radial 
coordinate; r t , thermal front; r c , char front; R, particle radius. 


called heat-balance integral method and was introduced by Good¬ 
man [20,21 ]. Hereafter, this phase of the biomass pyrolysis will be 
called heating of virgin biomass. 

When the temperature of the particle center starts to change, 
the second phase, called pre-pyrolysis heating, starts (Fig. 1(b)). 
The first two phases are different only in the mathematical for¬ 
mulation. The pyrolysis phase starts when the surface reaches the 
pyrolysis temperature. During this phase (Fig. 1(c)), the char front 
moves from the surface towards the center of the particle (shrink¬ 
ing core model). A further assumption is that the decomposition 
occurs instantaneously in the thin layer where the temperature 
equals the pyrolysis temperature. This allows to make use of an 
integral method, which transforms the partial differential equa¬ 
tions into ordinary differential equations. Haseli et al. [5] remark 
that, although the choice of a fixed pyrolysis temperature may be 
controversial, according to other studies [22] the process can be 
well described by assigning a suitable value for this temperature, 
depending on the process conditions, in such a way that the over¬ 
all energy balance agrees well with experimental results. Haseli 
assumes that during pyrolysis the particle does not shrink, thus the 
volume integrals involved are performed over a constant volume 
of the particle. The more severe is the torrefaction pretreatment, 
the better this assumption holds for the small torrefied biomass 
particles considered in this paper. 

Finally, when the particle is completely converted, the char par¬ 
ticle is further heated until it reaches the gas temperature. This 
phase is called post-pyrolysis heating (Fig. 1(d)). 

If the particle diameter exceeds a certain value, depending on 
initial temperature, pyrolysis temperature and heat flux, it may 
happen that pyrolysis starts before the thermal front has reached 
the particle center (Fig. 2(b)). Therefore, Haseli et al. [5,6] distin¬ 
guish thermally thin and thermally thick particles. The phases of 
a thermally thin and a thermally thick particle are illustrated in 
Figs. 1 and 2, respectively. In the following, first the model for a 
thermally thin particle will be introduced and then the model for a 
thermally thick particle. 


Fig. 2. Stages of a thermally thick biomass particle undergoing pyrolysis: heating of 
virgin biomass (a), thermally thick (b), pyrolysis (c), and post-pyrolysis heating (d). 


2.1. f. Heating of virgin biomass 

The heat equation in spherical coordinates (assuming spherical 
symmetry) governing the temperature of the particle in the first 
stage is: 


dr 
3 1 


,2 9T\ 
dr 


(1) 


with r is the radial coordinate, t the time, T(r, t) the temper¬ 
ature, and oib the (constant) thermal diffusivity of the virgin 
biomass. In Haseli’s approach [5,6] this partial differential equation 
is transformed to an ordinary differential equation for the surface 
temperature of the particle by assuming a parabolic temperature 
profile: 

|(qf(r t )) = 60 a B <j (2) 

with /(r t ) = 6 R 2 - qRr t + r t 2 + rf /R + rf/R 2 , r t the thermal front 
position, R the radius of the particle, and q the external heat flux: 

q = h(T g -T s ) + cre(T^-T s 4 ), (3) 


written as the sum of convective and radiative heat flux. The 
convective heat flux is proportional to the relative temperature 
(T g - T s ) between the surrounding gas and the particle surface. The 
heat-transfer coefficient h is computed using a correlation from 
literature that is applicable to spherical particles in the range of 
Reynolds (Re) and Prandtl (Pr) numbers relevant in our simulations. 
For forced convection around a sphere, the heat-transfer correla¬ 
tion is chosen as [23]: 

™ = i+0.3Rep /2 Pr'/ 3 , (4) 

where k g is the gas thermal conductivity, Re p = | v, - u(Xj, t) 12R/v 
is the Reynolds number based on the particle diameter and the 
relative velocity between the particle and the carrier gas at the 
particle position, and v is the kinematic viscosity of the gas. 

The radiative flux takes into account the heat flux received by 
the particle from the walls, whose temperature is indicated by T w 
and is constant. The Stefan-Boltzmann constant and the emissivity 
are indicated by a and €, respectively. 
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The thermal front position r t is expressed in terms of T s by eval¬ 
uating the parabolic temperature profile inside the particle at r=R: 
2ka 

r t = R- Z * B {T s - T 0 ) (5) 

with k B the thermal conductivity of the virgin biomass and To the 
initial particle temperature. Eq. (2) is the governing equation in the 
first phase. The unknown in (2) is T s . Eq. (2) can be written in the 
final form: 



with dTgldt the time derivative of the gas temperature at the parti¬ 
cle location. 

For times close to the initial time, the solution of this equation is 
approximately proportional to Vt. Therefore, explicit time integra¬ 
tion methods are not capable of solving this equation accurately. 
In order to circumvent this problem, we analyze the solution for 
small values of t in more detail. For small values of t, T s will be close 
to Tq and q can be assumed constant. A Taylor expansion of/[r t ) for 
T s close to T 0 shows that/(r t ) is proportional to ( T s - T 0 ) 2 . It follows 
that (T s - T 0 )~Vt. Therefore, we solve a differential equation for 
(T s — T 0 ) 2 , which equals: 

a . 2Cr s -r 0 )(6o aB q(r s -r„)-/^f 

f Ws+i-k (siM + d) 


2.1.2. Pre-pyrolysis 

The difference with the previous phase is the absence of the 
thermal front. Following Haseli et al. [5] the governing equation in 
this phase is: 


d_ 

di 


1 R 2 q 


(8) 


co = pci pb, with C vo i and C c the volatiles and char specific heat capac¬ 
ities, respectively. This term is not considered by Haseli in [5,6] but 
it has been taken into account in [25], 

This problem consists of two heat diffusion problems similar 
to the pre-pyrolysis phase, connected to each other at the inter¬ 
face r c . Because of the moving char front and the co-existence of 
these two regions, in the pyrolysis phase there are three governing 
ordinary differential equations, whereas in the other phases there 
is only one. The first governing equation is obtained from energy 
conservation in the char region and reads as follows: 

^/(r c , q, v c , 0i) = 60a c {R 2 q + k B r 2 fa -r 2 v c Ah p p B (13 j 

+(1 - a>)p B VcC vo ir 2 (T s - T p )], 

where /(r c , q, v c , 0i) = 9R 4 a + 6R 4 q - 11 R 3 ar c - 9R 3 qr c + R 2 qr 2 - 
R 2 ar 2 + Rqr 3 - Rar 3 + 4arf + qrf, a = v c p B Ah p - ks0i, with v c the 
velocity of the char front, and Ah p the specific heat of pyrolysis. 

Similarly, the energy conservation in the virgin biomass region 
yields the second governing equation in the pyrolysis phase: 

Tt(js (t) ^) = - 0lBr ^- (14) 


Finally, the third governing equation is the definition of the char 
front velocity: 



The unknowns in the three governing equations are T s , v c , and 
0i. The governing equations is written in matrix form. This set of 
ordinary differential equations is numerically integrated in time 
[24], 

2.1.4. Post-pyrolysis heating 

This phase is similar to the pre-pyrolysis phase with biomass 
properties replaced by char properties. Hence, the governing equa¬ 
tion is (9) with a B and k B replaced by a c and k c , respectively. 


which is solved for T s by numerically integrating the following 
equation: 



A detailed derivation of these equations is given in [24], 


2.1.5. Thermally thick particle 

The temperature profile inside a thermally thick particle con¬ 
sists of three regions (see Fig. 2(b)): 

T(r) = T 0 , for 0 < r < r t (16) 

T(r) = 0 2 (r-r c ) 2 +0i(r c -r) + 0o, for r t <r<r c (17) 

r(r) = 02(r-rc) 2 + 0i(r-r c ) + 0o, for r c <r<R, (18) 


2.1.3. Pyrolysis 

In the pyrolysis phase, the particle is divided into two regions 
separated by the interface at which pyrolysis occurs instanta¬ 
neously. The position of the interface r c is identified by the char 
front moving from the surface towards the center of the particle. In 
the model the temperature profiles in each region are assumed to 
be quadratic functions: 

T(r) = 0 2 (r- r c ) 2 + 0i(r c -r) + 0o, for 0 <r<r c (10) 

T(r) = 0 2 (r-r c ) 2 + 0 1 (r-r c ) + 0 o , for r c <r<R. (11) 


The heat equations governing this problem ; 
0<r<r c and: 


1 9 ( 
= ac r2Tr[ r 


,3T\ , rh 3T . 

S . f°rr t <r<K 


with etc the thermal diffusivity of the char. The second term on 
the right hand side of (12) represents heat transfer between the 
volatiles released at the char front and the char during their 
transport towards the particle surface. In (12) a V oi = Cvoi/PcQ and 


whose coefficients have to be determined using the boundary 
conditions. The energy conservation in the virgin biomass region 
r t < r < r c leads to the first governing equation in this phase: 

iKr t ,r c )=-60a B -^—, (19) 

at r c - r t 

with h(r t ,r c ) = 4r 3 + 3 r 2 r t + 2 r c r 2 + rf. The energy conservation in 
the char region r c <r<R leads to the second governing equation: 

^ g(tc , v c , r t ) = 60a c |r 2 9 - r 2 ( p B v v Ah p + 2k B T ^ _^° ) j (20) 

with g(r c , v c , r t ) = q(6R 4 - 9R 3 r c + R 2 r 2 + Rr 3 + r 4 ) + (2 k B (T p - 
ToVCtc ~ r t ) + p B v c Ah p )(9R 4 - 1 lR 3 r c - R 2 r 2 - Rr 3 + r c 4 ) + (1 - 
co)p B v c C vo ir 2 (T s — T p ). The third governing equation is (15), as for a 
thermally thin particle. 

2.2. Stiffness of the equations in the pyrolysis phase 

In a DNS with two-way coupling of mass, momentum and 
energy, it is essential that the properties of both phases (carrier 
and dispersed) are computed at the same time level. However, the 
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Fig. 3. Scheme of the routes of a biomass particle undergoing pyrolysis - in the flow diagram the criterion which is satisfied first dictates the route followed by the particle. 


system of equations for the particles in the pyrolysis phase can¬ 
not be solved with the same integration method as the one used 
for the Navier-Stokes equations of the gas. In the early and late 
phases of pyrolysis, the particle equations are stiff, for both the 
thermally thick and the thermally thin particle models. The stiffness 
problem can be circumvented by analysis of the problem for very 
small char and virgin biomass regions (see 3.3 in [24]). By applying 
approximations based on this analysis in the early and late phases 
of pyrolysis, the pyrolysis phase is divided into three sub-phases: 
early stiff pyrolysis sub-phase, regular pyrolysis sub-phase, and late 
stiff pyrolysis sub-phase. The equations of the second sub-phase are 
the same as presented before, for the thin and thick particle mod¬ 
els. The equations in the other two sub-phases can be derived in 
the same way as before by assuming a linear temperature profile in 
the char region in the early pyrolysis sub-phase and constant tem¬ 
perature profile in the virgin biomass region in the late pyrolysis 
sub-phase. 


2.3. Biomass particle routes undergoing pyrolysis 

We can distinguish two possible routes of a thermally thick 
particle undergoing pyrolysis, which, together with the route of 
a thermally thin particle, form a total of three possible routes for a 
pyrolyzing biomass particle. 

In order to have a better picture of the various phases, the routes 
of a biomass particle undergoing pyrolysis are schematically shown 
in Fig. 3. The scheme also shows the criteria used to switch between 
the phases. By e R a small value close to zero is indicated. Initially, the 
virgin biomass is heated. Afterwards, depending on the particle size 
and the heat flux applied, the particle may follow three different 
routes: 

1 A thin particle will move to the thin particle pre-pyrolysis phase 
when r t < Cr. Secondly, when T s = T p to the early stiff phase of pyrol¬ 
ysis and next, to the pyrolysis phase if r c < OR, with a typical value 
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of 6 = 0.9. Then, if <p\ < € ( j >t the equations of the late stiff phase of 
pyrolysis are solved. And finally, the particle moves to the post¬ 
pyrolysis or char heating phase. 

2 A thick particle will move to the thick particle early stiff phase 
when T s = Tp. Then, if r t < before r c <6R, the thick particle 
moves to the early stiff phase of pyrolysis and follows the same 
route as a thermally thin particle. 

3 Otherwise, the thermally thick particle moves to the thick par¬ 
ticle pyrolysis with moving thermal front and when r t < Er to the 
pyrolysis phase. Next, it follows the same route as a thin particle. 

In the following, the application of this biomass pyrolysis model 

in DNS is presented. 

2.4. Particle equation of motion 


We adopt the formulations for drag force on a particle and heat 
transfer between particle and gas for spherical particles. The same 
assumption is used by Papadikis et al. [26] and Ma et al. [27], 
who introduced a shape factor. Incorporation of the effects of a 
non-spherical particle shape on the correlations for drag force and 
convective heat transfer is a topic of future research. 

Since particles are small and have a much higher mass density 
than the gas, the drag force and the gravity are the dominant forces 
exerted on a particle [28], We do not take gravity into account since 
settling of particles is not a focal point in this study. Therefore, 
Newton’s law applied to a particle can be written as: 


d(m,v,) 

dt 


. (l + 0.15Rep 687 ) dnq 


(21) 


Here, m,- is the mass of particle i, v, its velocity and t p the par¬ 
ticle relaxation time given by r p = p p d p /(18/r, g ). Moreover, d p is 
the particle diameter, p p is the mass density of the particle and 
pig the dynamic viscosity of the gas. The location of the particle is 
denoted by x,- and u(x„ t ) is the velocity of the gas at x,. The two 
terms on the right-hand side are the drag force, where the stan¬ 
dard Schiller-Naumann drag correlation valid for particle Reynolds 
numbers between 0 and 1000 is adopted [29], and the change in 
momentum due to phase change. Eq. (21) can also be written in 
terms of accelerations, yielding the following particle equation of 
motion: 


dVj 

~dt 


' (l + 0.15Re p 687 ) 


(22) 


The Lagrangian particle tracking requires integration of the trajec¬ 
tory equation: 


dXj(f) 

dt 


(23) 


in which the particle velocity v, is computed from (22). 


3. Gas model 


The computational gas model in our simulations is composed 
of the compressible Navier-Stokes equation for momentum trans¬ 
port, coupled to the continuity equation, and the total energy 
equation. The dispersed phase (particles) is coupled to the carrier 
phase (gas) by means of source terms reflecting the presence of 
particles. The coupling between the phases can be classified into 
three types: (1) “one-way coupling” if only the flow influences the 
particles: (2) “two-way coupling” if particles have also a noticeable 
feedback on the flow; (3) “four-way coupling” when the particle 
concentration is so large that direct particle-particle interaction is 
important as well [10], We adopt two-way coupling, which involves 
momentum, mass and heat exchange between gas and particles. In 
this formulation, we assume the same properties for the volatiles 


and the nitrogen. Therefore, the mass lost by the particles is directly 
added to the gas, which is considered to be a single component gas, 
consistent with the assumption of low particle concentrations, and 
hence relatively low amounts of volatiles released compared to the 
total system. A model for the volatiles composition is not required 
in this work since only pyrolysis is considered. 

The source terms in the gas equations resulting from two-way 
coupling satisfy the requirement that they only transfer mass, 
momentum and internal energy from one phase to the other. They 
can be written as: 


/ v L2-dV--£s 


where the sum is taken over all particles within volume V. The first 
three components are the momentum exchange between the two 
phases, expressed by (21 ). The fourth component is the mass given 
by the particle to the gas during pyrolysis in the form of volatiles. 
The mass rate of change or mass loss is: 

^ = 4tt{pb - PcYcVc (25) 

In our model, the volatiles added to the gas phase cause an increase 
in the gas mass density and, consequently, in the pressure, since the 
volume of our domain is constant and periodic boundary conditions 
are applied in the streamwise and spanwise directions. For the sim¬ 
ulations with the highest mass load of biomass considered here, the 
resulting increase in pressure is almost 50%. The fifth component 
is the two-way coupling term of the energy, which is the sum of 
three contributions: 

te,.M,n-T.nc, J T r d ^-t (lm.lv/). 


The first term is the convective heat transfer between gas and 
particle, which is proportional to the temperature difference, the 
particle surface area A, and the convective heat transfer coefficient 
as defined by Bird et al. [23], The second term is the energy asso¬ 
ciated with the volatiles given to the gas by the particle during 
pyrolysis. The volatiles are formed at the char front where the tem¬ 
perature is T p . During their traveling to the particle surface, their 
increase in temperature is negligible. Therefore, to save computa¬ 
tional time, the temperature of the released volatiles is set to T p . 
In this work, the specific heat capacity of the volatiles and gas is 
taken equal. The third term is the change in kinetic energy which 
is expressed by the drag force and the mass loss of the particle: 


d_ 

dt 


dt 


in which dv,/dt is given by (22). 


4. Geometry and numerical method 

The computational domain is a channel with size 4 tzH in stream- 
wise direction and 2nH in spanwise direction, where H=0.2m 
is half the channel height (see Fig. 4). For the streamwise, wall- 
normal, and spanwise directions the notation x, y, and z is used, 
respectively. There are 128 3 control volumes in total. The flow 
solver is described in detail in [30], The spatial discretization is 
based on a finite volume method with second-order accuracy. 

Due to the coupling terms, it is important that the particle and 
the gas quantities are evaluated at the same time level. This is 
obtained using a low-storage second-order accurate Runge-Kutta 
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Fig. 4. Geometry of the channel. 

method for both gas and particle equations. The time step At is 
chosen small enough to ensure a stable and accurate solution of 
the equations. 

In the finite volume method for the gas, the two-way coupling 
terms (24) account for all particles present in the cell. The quanti¬ 
ties in (24) are already computed in the particle equations, which 
simplifies the implementation of these terms. Details on how to 
numerically treat the coupling terms are described in [30,24], 

5. Set up of the simulations 

Simulations are performed at a frictional Reynolds number 
approximately equal to Re r = 150, based on friction velocity u T = 
yjzw/Pg and half the channel height, where r w = ft g ^ | y= ±H is 
the wall shear stress and (u) the mean streamwise velocity com¬ 
ponent, averaged over the two homogeneous directions and time. 
The grid resolution is fine enough for capturing the main elements 
of the flow characteristics without employing a turbulence model. 
The initial volume fraction varies from simulation to simulation but 
it is always smaller than 0( 10 -3 ), so that particle-particle collisions 
can be ignored but two-way coupling might be relevant [10], In all 
simulations shown here the particle size is small compared to the 
smallest scales of turbulence, quantified by the Kolmogorov length, 
throughout the channel. Therefore, the point-particle approach 
employed in this work, is justified. 

In Table 1 the properties of the virgin biomass and the char are 
shown. The virgin biomass density is about 3.4 times larger than 
the char density. Therefore, a biomass particle loses most of its 
mass during conversion, releasing large amounts of volatiles. The 
constants in the radiative heat flux are a = 5.67 x 10 -8 W/(m 2 I< 4 ) 
and e = 0.9. The value of the specific heat of pyrolysis is 
Ah p = —270 kj/kg. 

The flow is initialized by a turbulent velocity field in the statisti¬ 
cally steady state obtained from a simulation without particles. The 
initial temperature of the gas is approximately constant at 1400K, 
whereas the initial temperature of the particles equals 300 K and is 
constant within the particle. The wall temperature is kept constant 
at 1400 K [31], Initially, particles are randomly and uniformly dis¬ 
tributed and the particle velocity is set equal to the gas velocity at 
the particle location. In the initial transient phase the temperature 
changes because of the heat flux q which transfers energy from 
the gas and from the walls to the colder particles by convection 


Properties of the virgin biomass and the biomass char used in the simulations. 

Mass density Thermal conductivity Specific heat 

(kg/m 3 )(W/(mK))capacity (J/(kgK)) 
Biomass 650 0.25 2500 

Char 190 0.1 1100 


and radiation. Using these settings we simulate an experiment of 
biomass particles undergoing pyrolysis in a turbulent channel flow 
with initial isothermal temperature and constant temperature of 
the walls. 


6. Results 

The conversion time of pyrolysis is the time required by the 
particles to be converted from biomass to char. It includes the time 
needed to heat the particles from their initial temperature to the 
pyrolysis temperature. In the results shown in this section, the 
conversion time is evaluated at 99% conversion of virgin biomass 
matter. Results evaluated at 95% of conversion are qualitatively the 
same. 

In Sections 6.2 and 6.3, we investigate the effect of varying vol¬ 
ume fraction and particle size on the conversion time. In Section 
6.4, we compare the results of the simulations with one-way and 
two-coupling for a realistic biomass particle size distribution. 

6.1. Effects of two-way coupling 

The effect of two-way coupling can be observed by studying the 
time evolution of the temperature profiles of nitrogen and particles. 
In Fig. 7 the mean nitrogen temperature and mean particle surface 
temperature are shown as functions of time. The volume fraction 0 
is 5 x 10 -7 . The particle diameter is 0.4 mm. Initially, the gas is much 
warmer than the particles, so that the particles receive energy from 
the gas and from the walls. The volume fraction is very small and 
the gas temperature is only slightly affected by the presence of the 
colder particles. The particle temperature rises until it reaches the 
gas temperature. Afterwards, the particle temperature becomes a 
little higher than the gas temperature due to the radiative heat flux 
from the walls, whose temperature is assumed constant (at 1400 K). 
From these results we conclude that the number of particles in this 
simulation is too small to appreciate two-coupling effects. Thus, for 
computing the conversion time of the biomass particle, a constant 
gas temperature could be a good approximation without requiring 
the employment of a numerical solver for the gas phase. 

A different scenario is presented in Fig. 8. The nitrogen and par¬ 
ticle surface mean temperatures are averaged in the homogeneous 
directions and plotted as functions of the wall-normal coordinate 
at different instants of time. In this simulation the diameter is the 
same as above but the volume fraction 0 is 2.55 x 10 -4 and two- 
way coupling effects cannot be neglected. In fact, the gas starts 
to cool down, in the center of the channel more than near the 
walls (the walls are kept at constant temperature). After a cer¬ 
tain time (around t= 0.17 s), particles become warmer than the gas 
because they are continuously heated by the radiative flux from 
the walls. Therefore, the convective heat exchange changes direc¬ 
tion and moves from the particles to the gas. This is due to the lack 
of a heat source for the gas apart from the heat diffused from the 
hot walls. Indeed, nitrogen does not allow any reaction to take place 
except for pyrolysis. When the particles are completely converted 
to char, the char heating continues until the particle temperature 
equals the wall temperature. 

In Fig. 5, isolines of the gas temperature in a vertical plane 
of the channel are shown during the initial heating phase of the 
particles. In the two-way coupling formulation (above), the par¬ 
ticles can be identified at the colder spots with lighter colors. In 
the one-way coupling formulation the particles do not affect the 
temperature field (figure below). Isolines of temperature during 
pyrolysis in Fig. 6 show the decrease in gas temperature due to the 
particles in the two-way coupling formulation, whereas in the one¬ 
way coupling formulation the temperature field is not significantly 
changed. 
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Fig. 5. Heating phase: instantaneous gas temperature in Kelvin in a vertical plane of the channel for the simulation with 0 = 2.55 x 10 4 and d p = 0.4 mm; above: two-way 
coupling formulation; below: one-way coupling formulation. 


We performed simulations with different numbers of particles 
in order to analyze the dependence of the pyrolysis conversion time 
on the particle volume fraction. 

6.2. Dependence on volume fraction 

The diameter of all particles was 0.4 mm. The results in Fig. 9 
show that in case two-way coupling is applied the conversion time 
increases with increasing volume fraction and it is constant for 
one-way coupling. The two-way particle-gas interaction affects the 
gas temperature because of the convective heat exchange: more 
particles extract more heat from the gas compared to simulations 
with fewer particles. This results in a slower particle heating and 
increases the pyrolysis conversion time. In the one-way coupling 
simulations, the temperature is not affected by the presence of par¬ 
ticles. In terms of heat exchange, the comparison between one-way 
and two-way coupling in Fig. 9 shows that the particle feedback to 
the gas becomes important for accurately predicting the conversion 


time when 0 > 10 -5 . The conversion time for 0 = 10 -3 is twice the 
conversion with one-way coupling, whereas for 0< 10 -4 the two- 
way coupling only leads to negligibly small differences. The range of 
volume fractions in which two-way coupling for the energy is rele¬ 
vant is equal to the range in which two-way coupling of momentum 
is important [10], According to Elghobashi [10], for 0> 10 -3 four¬ 
way coupling of momentum needs to be considered. We limited 
our investigation to two-way coupling so we cannot say whether 
or not the four-way coupling for the energy (radiation between the 
particles) is also important if 0 > 10 -3 . 

6.3. Dependence on particle size 

In order to investigate the effect of particle size on the conver¬ 
sion time, we performed simulations at constant volume fraction 
and varying particle diameter. Fig. 10 shows the conversion time 
tconv dependence on the particle diameter for three different vol¬ 
ume fractions: low, medium and high. For each volume fraction, 
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Fig. 7. Mean gas temperature (-) and particle surface temperature (—) varying in 
time for a low volume fraction (0 = 5 xlO -7 ). 


tconv increases with increasing diameter. For the same diameter, 
the conversion time t CO nv at a higher 0 increases more compared 
to a lower 0 or one-way coupling. The effect of two-way coupling 
is more clear using a logarithmic scale in which it is visible that, 
depending on the volume fraction, t CO nv is proportional to a cer¬ 
tain power of dp, whereas for the one-way coupling formulation 
(or very low 0) t CO nv ocR 5 / 3 . The figure also shows that the relative 
increase of t CO nv is larger for small particles. This means that the 
two-way coupling effect on the gas is reduced if, at same volume 
fraction, the particles are bigger. This can be explained by the larger 
total surface area of smaller particles compared to bigger particles 
with the same volume fraction and by the fact that the heat transfer 
coefficient is inversely proportional to the diameter of the particle 
(Eq. (4)). 

According to Fig. 8, the delay in the conversion time due to 
two-way coupling strongly depends on 0 but seems to be approx¬ 
imately independent of d p . These results suggests that, in first 
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approximation, the delay in the conversion time for a given diame¬ 
ter and volume fraction can be estimated from the conversion time 
computed with a one-way coupling formulation, which is much 
less computationally expensive. 

6.4. Particle size distribution 

In practical applications of co-firing of biomass and coal, the 
probability density function of biomass particle diameter follows a 
Rosin-Rammler distribution [32], This function P(d p ) is defined as 
follows: 

P{dp)= i(i) e ~ (dp/xf (26) 


with k and X the shape and scale factors, respectively. 

Here, we consider the size distribution of torrefied beech wood 
‘B-280-5’ from Repellin et al. [18], This distribution refers to beech 
chips torrefied for 5 min at 280 °C. The Rosin-Rammler distribution 
which fits these data is shown in Fig. 11 , with fc = 1.17, A. = 0.18 and 
the diameter in millimeter. The size of torrefied biomass is gen¬ 
erally smaller than the size of raw biomass. How small torrefied 
biomass particles depends on torrefaction conditions and inten¬ 
sity of milling [18], After the torrefaction, the biomass is lighter 
and easier to be ground. This yields smaller biomass particles after 
milling. 

The presence of very small particles requires a very small 
time step because of the stability restriction associated with the 
explicit time-integration method used here. Including very small 
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particles would lead to very long computational times. Therefore, 
we bounded the size distribution from below, in such a way that 
very small particles are disregarded without noticeably influencing 
the results [24], This saves computing time. 

In order to study the importance of gas-particle interaction on 
the pyrolysis conversion time, we compare the particle mass histo¬ 
ries of the two-way coupling simulation and a simulation without 
two-way coupling. The comparison in Fig. 12 suggests that two-way 
coupling effects are of minor importance due to the low particle vol¬ 
ume fraction used in these simulations. The delay is approximately 



Rosin-Rammler distribution (—) and experimental data (•) [18], 



t (s) 

Fig. 12. Normalized mass history of the size distribution of torrefied beech wood 
‘B-280-5’ from Repellin et al. [ 18 ] normalized by the initial mass at volume fraction 
= 5.23 x 10 -5 ; one-way coupling (—); two-way coupling (—). 



Fig. 13. Normalized mass history of the size distribution of torrefied beech wood 
‘B-280-5’ from Repellin et al. [18] normalized by the initial mass at volume fraction 
= 1.0 x 10- 3 ; one-way coupling (—); two-way coupling ( ). 

8% compared to the one-way coupling simulation. Indeed, the vol¬ 
ume fraction is 5.23 x 10 -5 , equivalent to 1000 kg of biomass for 
7419 kg of air, corresponding to 20% of air excess with respect to 
the stoichiometric particles-gas mass ratio (1000/6360) [31 ]. 

However, in our model combustion is not modeled. Therefore, 
the amount of biomass can be increased. Moreover, also with com¬ 
bustion in Oxyfuel conditions, where the oxygen concentration is 
very high, the biomass ratio can be higher than the stoichiometric 
value. Increasing 0 by a factor of 20 (0 = 1.044 x 10 -3 ) shows that 
the particle feedback to the gas delays the conversion time with 
more than 50% with respect to the same simulation with one-way 
coupling (Fig. 13). 

7. Conclusion 

In this paper we have shown how to extend Haseli’s pyrol¬ 
ysis model for biomass to an Euler-Lagrange formulation which 
makes the model suitable for implementation in a DNS code. Our 
formulation offers the advantage to directly simulate the interac¬ 
tion between gas and particles and to accurately predict particle 
dispersion. These two characteristics are very important in simu¬ 
lating biomass pyrolysis that may strongly depend on gas-particle 
interaction and particle dispersion, as we have shown in the 
results. In particular, the heat and mass exchange between par¬ 
ticles and gas influences the biomass pyrolysis conversion. The 
conversion time increases with increasing volume fraction and 
particle size. The particle feedback to the gas is only important 
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if 0>1O -5 . The conversion time delay is proportional to a certain 
power of the particle diameter depending on the volume fraction. In 
first approximation, the delay can be estimated from the one-way 
coupling result, thus saving computational power. Our formula¬ 
tion allows the simulation of a realistic biomass size distribution 
for predicting the biomass loss history. Theoretically, any size dis¬ 
tribution could be simulated and, based on the volume fraction, the 
particle-gas interactions are also taken into account. For the distri¬ 
bution of torrefied biomass considered in this work, it was found 
that the very small particles have a minor effect on the conversion 
time. Therefore, particles with a diameter smaller than 0.2 mm can 
be ignored. 
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